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Abstract 

We introduce a general distributional framework that results in a unifying description and character- 
ization of a rich variety of continuous-time stochastic processes. The cornerstone of our approach is an 
innovation model that is driven by some generalized white noise process, which may be Gaussian or not 
(e.g., Laplace, impulsive Poisson or alpha stable). This allows for a conceptual decouphng between the 
correlation properties of the process, which are imposed by the whitening operator L, and its sparsity 
pattem which is determined by the type of noise excitation. The latter is fully specified by a Levy 
measure. We show that the range of admissible innovation behavior varies between the purely Gaussian 
and super-sparse extremes. We prove that the corresponding generalized stochastic processes are well- 
defined mathematically provided that the (adjoint) inverse of the whitening operator satisfies some Lp 
bound for p > 1. We present a novel operator-based method that yields an explicit characterization of all 
Levy-driven processes that are solutions of constant-coefficient stochastic differential equations. When 
the underlying system is stable, we recover the family of stationary CARMA processes, including the 
Gaussian ones. The approach remains valid when the system is unstable and leads to the identification 
of potentially useful generalizations of the Levy processes, which are sparse and non-stationary. Finally, 
we show how we can apply finite difference operators to obtain a stationary characterization of these 
processes that is maximally decoupled and stable, irrespective of the location of the poles in the complex 
plane. 
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I. Introduction 

In recent years, the research focus in signal processing has shifted away from the classical linear 
paradigm, which is intimately linked with the theory of stationary Gaussian processes |[T|, Q. Instead of 
considering Fourier transforms and performing quadratic optimization, researchers are presently favoring 
wavelet- like representations and have adopted sparsity as design paradigm ||3|-||7]l. The property that a 
signal admits a sparse expansion can be exploited elegantly for compressive sensing, which is presently 
a very active area of research (cf. special issue of the Proceedings of the IEEE |[8|, Q). The concept is 
equally helpful for solving inverse problems and has resulted in significant algorithmic advances for the 



efficient resolution of large scale £i-norm minimization problems |10|-|12|. 

The current formulations of compressed sensing and sparse signal recovery are fundamentally deter- 
ministic. By drawing on the analogy with the classical theory of signal processing, it is likely that further 
progress may be achieved by adopting a statistical (or estimation theoretic) point of view. This stands as 
our primary motivation for the investigation of the present class of continuous-time stochastic processes, 
the greater part of which is sparse by construction. These processes are specified as a superset of the 
Gaussian ones, which is essential for maintaining backward compatibility with traditional statistical signal 
processing. 

The inspiration for this work is provided by the innovation approach to system modeling — a standard 
technique in statistics and control theory that is well developed in the discrete setting and often favored by 
engineers. Innovation models are also used in signal processing for the investigation of continuous-time 
stationary Gaussian stochastic processes ||T|, p3|. Non-Gaussian variants of such models are easy to set 



up in the discrete world, but they do result in harder identification problems |14|-|16|. An active topic 
of research is the determination of a proper noise input to simulate signals with prescribed marginal 



distributions |17|, 1 18|. By contrast, there is comparatively little work on continuous-domain innovations 
for the specification of non-Gaussian or/and non-stationary processes due to the inherent difficulty 
of rigorously defining non-Gaussian white noise in the continuous domain. The proper mathematical 
framework exists and was developed by the Russian school of mathematics in the 1960s [19], but has 
hardly been used by practitioners until now. This is mainly due to the widespread acceptance of stochastic 



integration (Ito calculus) in the advanced theory of stochastic processes pO|-p3|, which avoids the direct 
handling of white noise and tempered distributions. 

By following up on our initial work on the generation of piecewise-smooth signals from random 



streams of Dirac impulses (Poisson white noise) |24|, our present aim is to set the foundations of a 
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comprehensive theory of continuous-domain stochastic processes based on the simple, unifying principle 
of the filtering of special brands of (non-Gaussian) white noise. While the concept remains applicable 
in multiple dimensions, we focus on the time domain (1-D signals), and provide a systematic treatment 
of systems that are described by ordinary differential equations, including some novel twists for the 
non-stable scenarios, which opens the door to interesting generalizations. The primary contributions are: 
1) The extension of our prior innovation models to the broadest possible class of white noises beyond 
the Gaussian and impulsive Poisson categories: We show that each brand is uniquely specified 
by a Levy measure that conditions the degree of sparsity of the process. The Gaussian processes 
are the least sparse ones; the Poisson processes are intermediate with their level of sparsity being 



controlled by the rate parameter A |24|. The sparsest processes are the alpha-stable ones whose 



marginal distributions are heavy tailed with unbounded variance |22|, |25|. 

2) The systematic investigation of processes that are ruled by constant-coefficient SDEs together with 
the proposal of a generic operator-based method of solution: When the underlying system is stable, 
we recover the complete family of (non-Gaussian) continuous-time autoregressive moving average 
(CARMA) processes (see also the work of Brockwell for an equivalent state-space characterization 
that relies on stochastic integrals |26]). The further reaching aspect of our formulation is that the 
method remains applicable in the non-stable case and that it leads to some interesting generalizations 
of Levy processes, which are non-stationary. 

3) The generalization/extension of our previous stability and existence results (cf. pA\ Theorem 2], 



1 27 Theorem 1.3]) for our enlarged class of stochastic processes: In essence, we are replacing the 
basic L2-boundedness requirement that is central to the continuous-time Gaussian theory by a more 
robust Lp condition (cf. Theorem 3); the case p = 1 is required for the non-symmetric Poisson 
processes, while the range of values p G (0, 2) becomes appropriate for the alpha stable processes. 
The paper is organized as follows. In Section II, we briefly review the foundations of Gelfand's 
theory of generalized stochastic processes. In particular, we characterize the complete class of admissible 
continuous-time white noise processes and give some argumentation as to why the non-Gaussian brands 
are inherently sparse. We then introduce our general innovation model in Section III, while providing a 
novel operator-based method for the solution of SDE. In Section IV, we make use of Gelfand's formalism 
to fully characterize our extended class of (non-Gaussian) stochastic processes including the special cases 
of CARMA and A^th-order generalized Levy processes. We also introduce the notion of generalized 
increment process and some corresponding B-spline calculus, which allows for a common stationary 
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treatment of the two latter classes of processes, irrespective of any stability consideration. Finally, in 



Section V, we take a non-conventional look at the classical Levy processes |23|, |28|, |29|, which 



correspond to the simplest (unstable) instance of our extended formulation (single pole at the origin). 



For higher-order illustrations of sparse processes, we refer to our companion paper |30|, which is more 
specifically devoted to the study of the discrete-time implication of the theory. The notation, which is 
common to both papers, is summarized in |30J Table I]. 

II. Preliminaries 

The purpose of this section is to introduce the distributional formalism that is required for the proper 
definition of continuous-time white noise. We start with a brief summary of some required notions in 
functional analysis, which also serves us to set the notation. We then introduce the fundamental concept 
of characteristic functional which constitutes the foundation of Gelfand's theory of generalized stochastic 
processes. We proceed by giving the complete characterization of the possible types of continuous-domain 
white noises — not necessarily Gaussian — which will be used as universal input for our innovation models. 
We conclude the section by showing that the non-Gaussian brands of noises that are allowed by Gelfand's 
formulation are intrinsically sparse, a property that has not been emphasized before (to the best of our 
knowledge). 



A. Functional and distributional context 

The Lp-norm of a function f{t) is ||/||p = (/j^ dt) ^ for 1 < p < oo and ||/||oo = esssup^gR 

for p = +00 with the corresponding Lebesgue space being denoted by Lp. The concept is extendable for 
characterizing the rate of decay of functions. To that end, we introduce the weighted Lp ^ spaces with 

a G ]R+ 

Lp,a = {fit) G Lp : ||/||p,o < +00} 
where the a-weighted Lp-norm of / is defined as 

ll/IU = ll(i + ir)/Wllp- 

Hence, the statement / G Loo,« implies that f{t) is decaying at least like (l/|t|")ast tends to ±00; more 
precisely, that \ f{t)\ < almost everywhere. In particular, this allows us to infer that L^^ C Lp 

for any e > and p > 1. Another obvious inclusion is Lp q, C Lp for any a > ao- In the limit, we 
end up with the space of rapidly-decreasing functions TZ = {f{t) : ||/||oo,m < +00, Vm G Z+}, which 
is included in all the others. 
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We use ip{t) to denote a generic function in Schwartz's class S of rapidly-decaying and infinitely- 
differentiable test functions. Specifically, Schwartz's space is defined as: 

S = {^(t) G C7~ : ||D"99||oo,m < +00, Vm, n G Z+} , 

with the operator notation D*^ = and the convention that D*^ = Id (identity). 5 is a complete 
topological vector space. Its topological dual is the space of tempered distributions S'; a distribution 
€ 5' is a continuous linear functional on S that is characterized by a duality product rule ^(v?) = 
((/>,(/?) = f^(j){t)(p{t) dt with if ^ S where the right-hand side expression has a literal interpretation as 
an integral only when (l){t) is true function of t. The prototypical example of a tempered distribution 
is the Dirac distribution 5, which is defined as 5{lp) = {6,ip) = (p{0). In the sequel, we will drop the 
explicit dependence of the distribution on the generic test function ip £ S and simply write (p or even 
(with an abuse of notation). 

Let T be a continuou^ linear operator that maps S into itself (or eventually some enlarged topological 
space such as Lp). It is then possible to extend the action of T over S' (or an appropriate subset of 
it) based on the definition (Tc/;, = {(p,T*ip) if T* is the adjoint of T which maps to another test 
function T*ip ^ S continuously. An important example is the Fourier transform whose classical definition 
is T{f}{uj) = f{uj) = /^/(i)^ Since is a self-adjoint 5-continuous operator, it is extendable 

to S' based on the adjoint relation {Tcl),(p) = {(pjTip) for all ip £ S (generalized Fourier transform). 

A linear, shift-invariant (LSI) operator that is well-defined over S can always be written as a convolution 
product: 

Tl,sMt) = {h*ip)it)= [ h{T)^{t-T)dT 

JK 

where h{t) = TLsi5(t) is the impulse response of the system. The adjoint operator is the convolution 
with the time-reversed version of h: 

h'^it) = h{-t). 

The better-known categories of LSI operators are the BIBO-stable (bounded input, bounded output) filters, 
and the ordinary differential operators. While the latter are not BIBO-stable, they do work well with test 
functions. 

'An operator T is a continuous map from a separable topological vector space V into another one iff. ipk ^ (p in the topology 
of V implies that Tifk — >■ T<^ in the topology (or norm) of the second space. If the two spaces coincide, we say that T is 
V-continuous. 
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1) Lp-stable LSI operators: The BIBO-stable filters correspond to the case where h E Li, or, more 
generally, when h corresponds to a complex-valued Borel measure of bounded variation. The latter 
extension allows for discrete filters of the form hd{t) = J^nez d[n]5{t — n) with d[n] G ^i. We will refer 
to these filters as Lp-stable because they are bounded in all Lp-norms (by Young's inequality). Lp-stable 
convolution operators satisfy the properties of commutativity, associativity, and distributivity with respect 
to addition. 

2) S -continuous LSI operators: For an Lp-stable filter to yield a Schwartz function as output, it is 
necessary that its impulse response (continuous or discrete) be rapidly-decaying. In fact, the condition 
h eTZ (which is much stronger than integrability) ensures that the filter is iS-continuous. The nth-order 
derivative D" and its adjoint D"* = (— 1)"D" are in the same category. The nth-order weak derivative 
of the tempered distribution ^ is defined as D"(^((/?) = (D"^, tp) = {<p, D"'*ip) for any ip e S. The latter 
operator — or, by extension, any polynomial of distributional derivatives P/v(D) = Z]^=i "^nl-*" with 
constant coefficients a„ € C — maps S' into itself. The class of these differential operators enjoys the 
same properties as its classical counterpart: shift-invariance, commutativity, associativity and distributivity. 

B. Notion of generalized stochastic process 

The leading idea in distribution theory is that a generahzed function cj) is not defined through its point 
values <^{t),t G M, but rather through its scalar products 0((^) = (0, with all "test" functions (p e S. 
In an analogous fashion, Gelfand defines a generalized stochastic process s via the probability law of its 
scalar products with arbitrary test functions (p G S, rather than by considering the probability law of its 
pointwise samples {. . . , s{ti), s{t2), . ■ ■ , s{tN), ■ ■ ■}, as is customary in the conventional formulation. 

Let s be such a generalized process. We first observe that the scalar product xi = (5,991) with a 
given test function 991 is a conventional (scalar) random variable that is characterized by a probabil- 
ity density Pi{dxi); the latter is in one-to-one correspondence (via the Fourier transform) with the 
characteristic function pi{uJi) = S{e^'^^''^} = j^e^'^^''^ Pi{dxi) = S'{e^(^''^^'fi^')} where #{•} is tiie 
expectation operator. The same applies for the probabihty density Pi^2{dxidx2) associated with a pair 
of test functions (pi and ip2 which is the inverse Fourier transform of the 2-D characteristic function 
^i,2(^^i,'^2) = £'{e^^^''^'^'^'^'^"^'^^'>}, and so forth if one wants to specify higher-order dependencies. 

The foundation for the theory of generalized stochastic processes is that one can deduce the complete 
statistical information about the process from the knowledge of its characteristic form 

Z,(<^) = ,r{e^<^''^)} (1) 
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which is a continuous, positive-definite functional over S such that Zs{0) = 1. Since the variable ip in 
Zs{ip) is completely generic, it provides the equivalent of an infinite-dimensional generalization of the 
characteristic function. Indeed, any finite dimensional version can be recovered by direct substitution of 
if = uji^pi + • • • + ojn^pn in 2,s{^p) where the are fixed and where u = (wi, • • • , u^) takes the role 
of the A^-dimensional Fourier variable. In fact, Gelfand's theory rests upon the principle that specifying 
an admissible functional Zs{ip) is equivalent to defining the underlying generalized stochastic process 
(Bochner-Minlos theorem). The precise statement of this result, which relies upon the fundamental notion 
of positive-definiteness, is given in Appendix I. 

C. White noise processes 

We define white noise as a generalized random process that is stationary and whose measurements 
for non-overlapping test functions are independent. A remarkable aspect of the theory of generahzed 
stochastic processes is that it is possible to deduce the complete class of such noises based on functional 
considerations only p9| . To that end, Gelfand and Vilenkin consider the generic class of functionals of 
the form 

ZM = ^^vfyjj{v{t))dt^ (2) 

where f{u) is a continuous function of the scalar variable u G M. They first argue that this functional 
specifies an independent noise process iff. Zyj is positive-definite (as required by the Bochner-Minlos 
theorem) and Zw{lpi + (/?2) = Zyj{Lpi)Zyj{Lp2) whenever ipi and have non-overlapping support (i.e., 
• (/92(t) = 0). The latter is equivalent to the requirement that f{u) in ([2]) is such that /(O) = 0. They 
then go on to prove that the complete class of functionals of the form (|2]) with the required mathematical 
properties (continuity, positive-definitess and factorizability) is obtained when 

f{u) = bo+ jb[u - ^ + / [e^'^" - r(a)(l + jau)] R{da) 
^ Jm.\{o} 

where R is an other arbitrary positive measure on R such that /][j\^|o} niin(l, a^) R{da) < oo and where 
r(a) is a function such that r(a) — 1 has a third-order zero at a = and Jj^j^ot-'^ ~ ^('^)] R{da) + 6o = 
(to ensure that /(O) = 0). In doing so, they actually establish a one-to-one correspondence between the 
characteristic functional of an independent noise processes of the form ([2]) and the (classical) characteristic 



function e-^'^'^^ = S{e^^^} of an infinite divisible scalar random variable x |31|, |32|. An equivalent, 
more concise specification of the complete family of admissible functions f{u) is provided by the Levy- 
Khintchine formula 

f{u) = jb'lu - ^ + / [e^'^" - 1 - janl|,|<i(a)] V{da) (3) 



September 1, 2011 DRAFT 



8 



where l|(j|<i(a) is the indicator function that takes the value 1 if \a\ < 1 and zero otherwise. Here V is 
a Borel measure on M\{0} such that 



/ 

^R\{0} 



min(l, a^) V{da) < oo. (4) 
We will see that the density V{da) actually corresponds to the so-called Levy measure that enters the 



definition of Levy processes p3| , | |29| , |32|. To further our mathematical understanding of the Levy- 
Khintchine formula (jsjl, we note that e^°'^ — l—jaul\a\<:i{a) = —■^a'^u'^+0{\au\^) as a — )• 0. This ensures 
that the integral is convergent even when V{da) exhibits a singularity at the origin to the extent allowed 
by the admissibility condition (j4]). If the Levy measure is finite (i.e., V{da) < oo) or symmetrical 
iV{E) = V{—E) for any E C M), it is then also possible to use the equivalent, simplified form of Levy 
function 

,2 



f{u) = jhu - ^ + / (e^"" - 1) V{da) 



with h\ = 6" — l^^Y oy{do)- Concretely, this means that a particular brand of such independent noise 
process is completely characterized by a triplet (61,62,^(^0))- With this latter convention, the Levy 
functions and characteristic forms of the three primary types of white noise encountered in the signal 
processing Uterature are: 

1) Gaussian: 61 = 0, 62 = 1, ^ = 

/Causs ("w) = 1 

^«,(<^) = e-^ll'^ll^^. (5) 

2) Compound Poisson: h\ = 0, 62 = 0, V{dd) = A Pa(a) d« with J"jgPa(a) da = Pa(0) = 1, 

/Poisson(w; A, Pa) = / (c''"" " l) Vaifl) da, 



exp f A / / (e^'^'^(*) - 1) va{a) da At] . (6) 

3) Symmetric alpha-stable (SaS): 61 = 0, 62 = 0, V{da) = da with < a < 2 and = "'"^^"^^ 

a suitable normalization constant, 



Z^{ip) = e--^\\'^\\-^. (7) 



The latter follows from the fact that is the generalized Fourier transform of with the convention 



that a! = r(a + 1) where F is Euler's Gamma function [33J . 
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While none of these noises has a classical interpretation as a random function of t, we can at least 



provide an explicit description of the Poisson noise as a random sequence of Dirac impulses (cf. |24 
Theorem 1]) 

w'aW = '^akS{t - tk) 

k 

where the tk are random locations that are uniformly distributed over M with density A, and where the 
weights Ofc are i.i.d. random variables with probability density function (PDF) Pa{a). 

D. Gaussian versus sparse categorization 

To get a better understanding of the underlying class of objects, we propose to probe them through 
some localized analysis window ip, which will yield a conventional i.i.d. random variable x = {w,ip) with 
some probability density function (PDF) Pip{x). The most convenient choice is to pick the rectangular 
analysis window <p{t) = rect(t) = By using the fact that e^'^'^'^^^^i^) _ i = - 1 for 

t G [— ^, |], and zero otherwise, we find that the characteristic function of x is simply given by 

Prect(w) = Zw {u ■ rect(t)) = exp {f{uj)) , 

which corresponds to the generic (Levy-Khinchine) form associated with an infinitely-divisible distribution 



1 29 1, 1 32 1, |34|. The above result makes the mapping between generalized white noise processes and 
classical infinite-divisible (id) law^ explicit: The "canonical" id distribution of w, Pid{x) = Pvcct{x), is 
obtained by observing the noise through a rectangular window. Conversely, given the logarithm of the 
characteristic function of an id distribution, log (J^{pid(x)}(a;)) = /(w), we can specify a corresponding 
generalized white noise process w via the characteristic form Z^i^) by merely substituting the frequency 
variable co by the generic test function ip{t), adding an integration over M and taking the exponential as 
in 

We will now argue that this class of models allows for a range of behaviors that varies between the 
purely Gaussian and sparse extremes. In the context of Levy processes, these are often referred to as the 
diffusive and jump modes. To make our point, we consider two distinct scenarios. 

1) Finite variance case: We first assume that the second moment m2 = Jg^a^Vlda) of the Levy 
measure F in Q is finite. This allows us to rewrite the classical Levy-Khinchine representation as 

f{u) = jc,u - ^ + / [e^"" - 1 - jua] V{da) 
^ Jm.\{o} 

random variable x with PDF p{x) is said to be infinitely divisible (id) if for any n G N"*" there exist i.i.d. random variables 
xi, . . . ,x„ with PDF say Pn{x) such that x has the same distribution sls xi + ■ ■ ■ + x„. 
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where the Poisson part of the functional is now fully compensated. Indeed, we are guaranteed that the 
above integral is convergent because \e^"-^ — 1 — jua\ = 0{a^) as a — )• and je-'"" — 1 — jua\ = 0{\a\) 
as a — )• iboo. An interesting non-Poisson example of infinitely-divisible probability laws that falls into 
this category (with non-finite V) is the Laplace density with Levy triplet (0, 0, V{da) = da) and 
p{x) = ^e~l^'L This model is particularly relevant for sparse signal processing because it provides a tight 
connection between Levy processes and total variation regularization |24[ Section VI]. 

Now, if the Levy measure is finite V{da) = X < oo, the admissibility condition yields a V{da) < 
oo, which allows us to pull the bias correction out of the integral. The representation then simplifies to 

/(n) = jhu - ^ + A / [e^-- - 1] pa{a) da, 

with V(da) = Xpa{a) da and f^pa{a) da = 1. This implies that we can decompose x into the sum of 
two independent Gaussian and Compound-Poisson random variables. The variances of the Gaussian and 
Poisson components are cr^ = 62 and Ai?{a^}, respectively. The Poisson component is sparse because its 
probability density function exhibits a mass distribution e^'^6{x) at the origin, meaning that the chances 
for a continuous amplitude distribution of getting zero are overwhelmingly higher than any other value, 
especially for smaller values of A > 0. It is therefore justifiable to use < e-^ < 1 as our Poisson 
sparsity index. 

2) Infinite variance case: We now turn our attention to the case where the second moment of the 
Levy measure is unbounded, which we like to label as the "super-sparse" one. To substantiate this claim, 
we invoke the Ramachandran- Wolfe theorem which states that the pth moment with p G of 



an infinitely divisible distribution is finite iff. Jj^|>]^ \a\^ V{da) < 00 |35|, |36|. For p > 2, the latter is 
equivalent to Jj^ |a|^ V{da) < 00 because of the Levy admissibiUty condition. It follows that the cases 
that are not covered by the previous scenario (including the Gaussian -1- Poisson model) necessarily give 
rise to distributions whose moments of order p are unbounded for p >2. The prototypical representatives 
of such heavy tail distributions are the alpha-stable ones or, by extension, the broad family of infinite 
divisible probability laws that are in their domain of attraction. It has been shown that these distributions 
precisely fulfill the requirement for ip compressibility |37| , which is a stronger form of sparsity than the 
presence of a mass probability density at x = 0. 

III. Innovation approach to continuous-time stochastic processes 

Specifying a stochastic process through an innovation model (or an equivalent stochastic differential 
equation) is attractive conceptually, but it presupposes that we can provide an inverse operator (in the 
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form of an integral transform) that transforms the white noise back into the initial stochastic process. This 
is the reason why we will spend the greater part of our effort investigating suitable inverse operators. 

A. Stochastic differential equations 

Our aim is to define the generalized process with whitening operator L and noise parameters (61, 62, V) 
as the solution of the stochastic linear differential equation 

Ls = w, (8) 



where w is a white noise process, as described in subsection II-C This definition is obviously only usable 
if we can construct an inverse operator T = that solves this equation. For the cases where the inverse 
is not unique, we will need to select one preferential operator, which is equivalent to imposing specific 
boundary conditions. We are then able to formally express the stochastic process as a transformed version 
of a white noise 

s = L-^w. (9) 

The requirement for such a solution to be consistent with ([8]) is that the operator satisfies the right- 
inverse property LL^^ = I over the underlying class of tempered distributions. By using the adjoint 
relation {s,ip) = {L^^w,ip) = {w,L^^*ip), we can then transfer the action of the operator onto the test 
function inside the characteristic form and obtain a complete statistical characterization of the so-defined 
generalized stochastic process 

Zsiif) = Z^-.U^) = Z^(L-i», (10) 



where is given by (|2]) (or one of the specific forms in the list at the end of section II-C I and where 
we are implicitly assuming that the adjoint L~^* is mathematically well-defined over S. 

In order to fulfill the above mathematical requirements, it is usually easier to proceed backwards: one 
specifies an operator T that satisfies the left- inverse property: \/ip G S, XL* 99 = ip, and that is continuous 
(i.e., bounded in a proper metric) over the chosen class of test functions. One then characterizes the adjoint 
of T, which, for a given (j) e S, is such that 

Finally, one applies a standard limit argument to extend the action of T* = L^^ over the enlarged class 
of tempered distribution cj) £ S' based on the above adjoint relation, which yields the proper distributional 
definition of the right inverse of L in Q. 
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B. Inverse operators 

Before presenting our general method of solution, we need to identify a suitable set of elementary 
inverse operators that satisfy the required boundedness conditions. 

Our approach relies on the factorization of a differential operator into simple first-order components 
of the form (D — a-n^d) with a„ € C, which can then be treated separately. Three possible cases need 
to be considered. 

1) Causal-stable: Re(a„) < 0. This is the classical textbook hypothesis which leads to a causal- 
stable convolution system. It is well known from linear system theory that the causal Green function of 

(D — Unld) is 

where u{t) = l[o,+oo)(0 is the unit-step (or Heaviside) function. Clearly, Pa„{t) is absolutely integrable 
(and rapidly-decaying) iff. Re(a„) < 0. It follows that (D — a„Id)~^/ = pa„ * f with pa„ G TZ C Li. 
In particular, this implies that T = (D — a„Id)~^ specifies a continuous LSI operator on S. The same 
holds for T* = (D - Q„Id)-i*, which is defined as T* f = p\^ * f. 

2) Anti-causal stable: Re(a„) > 0. This case is usually excluded because the standard Green function 
Pa„{t) = n(t)e°"* grows exponentially, meaning that the system does not have a stable causal solution. 
Yet, it is possible to consider an alternative anti-causal Green function p'a^{t) = —P-a,S^) ~ /^"-(O ~ 
e""*, which is unique in the sense that it is the only Green functiorj^of (D — Onld) that is Lebesgue- 
integrable and, by the same token, the proper inverse Fourier transform of -^^^ for Re(Q„) > 0. In 
this way, we are able to specify an anti-causal inverse filter (D - a„Id)"V = Pa„ * / with p'^^^ £ U 
that is Lp-stable and 5-continuous. In the sequel, we will drop the ' superscript with the convention that 
Pa„{t) systematically refers to the unique Green function of (D — a„Id) that is in TZ (rapidly-decaying); 
it is either causal or anti-causal depending on the polarity of Re(a„)- 

3) Marginally stable: Re(a„) = or, equivalently, an = j^^o with uq G M. This third case, which is 
incompatible with the conventional formulation of stationary processes, is most interesting theoretically 
because it opens the door to important extensions such as Levy processes and fractals (fractional Brownian 
motion). Here, we will show that marginally-stable systems can be handled within our generalized 
framework as well, thanks to the introduction of appropriate inverse operators. 

^: pis a Green functions of (D — Q:„Id) iff. (D — Q„Id)~^p — 5; the complete set of solutions is given p{t) = pci„ (i) + Ce°"' 
which is the sum of the causal Green function pa„ (t) plus an arbitrary exponential component that is in the null space of the 
operator. 
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The first natural candidate for (D — ja;oId) ^ is the inverse filter whose frequency response is 

Pjujoi^) = -r-^ r + 7r5(w - Wo). 

It is a convolution operator whose time-domain definition is 

= (/5j<^o*V)(0 

= e-^'^"* / e-^'^''^(/^(r) dr (11) 

J — oa 

where Pjuig{t) = u{t)e^^°^; its adjoint is given by 

/+00 
e^'^^Xr) dr. (12) 

While I^^9?(t) and I*^^(/?(f) are both well-defined pointwise when ip £ Li, the problem is that these 
inverse filters are not BIBO stable since their impulse responses, Pjuoii) ^rid /'jajo(^)' ^i- 
particular, one can easily see that luj„(p (resp., I^^v?) with (/9 G 5 is generally not in Lp with 1 < p < +oo, 
unless ip{iOo) = (resp., ip{—ujQ) = 0). The conclusion is that I*^^ fails to be a bounded operator over 
the class of test functions S. 

This leads us to introduce some "corrected" version of the adjoint inverse operator I*^ 

Co.toV'W = C - 0i-^o)e-^^°'"6{- - to)} it) 

= ZMi) - '^(-^o)e-^''^'^*V)'.o(i - to), (13) 

where to G M is a fixed location parameter and where (^(— wo) = J^e^^"*Lp{t) dt is the complex 
sinusoidal moment associated with the frequency ujq. The idea is to correct for the lack of decay of 
1^0 v(*) as t — )• — oo by subtracting a properly weighted version of the impulse response of the operator. 
An equivalent Fourier-based formulation is provided by the formula at the bottom of Table I; the main 
difference with the corresponding expression for 1^1^'^ is the presence of a regularization term in the 
numerator that prevents the integrant from diverging at a; = wo- The next step is to identify the adjoint 
of I* ^ , which is achieved via the following inner-product manipulation 

{^Mao'i^) = {^M,<i>) - <^(-a;o)e-^'"»*° {^^p)u.S- - to)) (^Y linearity) 

^ V ' 

= iU^, ct)) - {e^^"\ (^) e-^''^"*" L„(/?(to) (using ^) 
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Since the above is equal to {li^g^to^j 'P) by definition, we obtain that 

i^aMt) = U^it)-e^''''^'-'"H^Mto). (14) 

Interestingly, this operator imposes the boundary condition l^g^taVi'^o) = via the substraction of 
a sinusoidal component that is in the null space of the operator (D — jujQld), which gives a direct 
interpretation of the location parameter to- Observe that expressions ( [T3] ) and ([14]) define linear operators, 
albeit not shift- invariant ones, in contrast with the classical inverse operators 1^^;^ and I*^^. 

For analysis purposes, it is convenient to relate the proposed inverse operators to the anti-derivatives 
corresponding to the case uoq = 0. To that end, we introduce the modulation operator 

which is a unitary map on L2 with the property that M~^^ = M^^^. 

Proposition 1: The inverse operators defined by ( [TT] ), (12i, (14i, and (13 1 satisfy the modulation 
relations 

Proof: These follow from the modulation property of the Fourier transform (i.e, F{yii^^,^p}{uj) = 
F{ip]{u) - u)q)) and the observations that = Pj^^{t) = M^^pQ{t) and = p)^S'^) = 

M_j^„/9^(t) with po{t) = u{t) (unit step). ■ 
The important functional property of I*^^ is that it essentially preserves decay and integrability, while 
lujg^ta fully retains signal differentiability. Unfortunately, it is not possible to have the two simultaneously 
unless luigip{to) and (p{—ujo) are both zero. 

Proposition 2: If / G Loo(i^,Wa) with a > 1, then there exists a constant Ct,, such that 

\iLtJ{t)\ < Ct,- 



1 + 

which implies that G i^oo,Q-i- 

Proof. Since modulation does not affect the decay properties of a function, we can invoke Proposition 
1 and concentrate on the investigation of the anti-derivative operator Iq j . Without loss of generality, we 
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can also pick to = and transfer the bound to any other finite value of to by adjusting the value of the 
constant Ctg. Specifically, for t < 0, we write this inverse operator as 

/■H-oo roo 

= / fir) dr - / fir) dr 
= - f fir) dr. 

This implies that 



/_JMd.|<||/IU„£^d.<(J^) 



1/ l + \t\ 

For t > 0, lQQ/(t) = /(r) dr so that the above upper bounds remain valid. ■ 
The interpretation of the above result is that the inverse operator I* ^ reduces inverse polynomial 
decay by one order. Proposition 2 actually implies that the operator will preserve the rapid decay of the 
Schwartz functions which are included in Loo,a for any a G It also guarantees that lZa,to^ belongs 
to Lp for any Schwartz function ip. However, I*^ will spoil the global smoothness properties of ip 
because it introduces a discontinuity at to, unless (^(— wq) is zero in which case the output remains in 
the Schwartz class. This allows us to state the following theorem which summarizes the higher-level part 
of those results for further reference. 



Theorem 1: The operator I* ^^ defined by (14 1 is a continuous linear map from S into TZ (the space 



of bounded functions with rapid decay). Its adjoint lu}o,to is given by (13 1 and has the property that 
It<;o,to9^(^o) = 0. Together, these operators satisfy the complementary left- and right-inverse relations 

i:o,io(D-i^oId)V = (/p 

(D - ju;oId)Lo,to93 = (f 

for all if € S. 

Having a tight control on the action of I*^ over S allows us to extend the right-inverse operator lujB,to 
to an appropriate subset of tempered distributions (j) G S' according to the rule {IujqM'P^ = Kjo U^)- 
Our complete set of inverse operators is summarized in Table I together with their equivalent Fourier-based 
definitions which are also interpretable in the generalized sense of distributions. 

C. Solution of generic stochastic differential equation 
We now have all the elements to solve the generic stochastic linear differential equation 

N M 

J]a„D-s= ^6^D"u; (15) 

n=l m=l 



September 1, 2011 



DRAFT 



16 



TABLE I 

First-order differential operators and their inverses 



Properties of inverse operator 



Standard case: a„ G C, Re(Q,i) 7^ 
(D - a„Id) 

(D - Q„Id)* 

Critical case: a„ — jujoj^o G K 
(D - jujold) 



(D*-a„Id)-7(f)= / /H 

Jr \-juj-ar, 



2tt 



duu 



J(t) = / /(c.) f — ^ 
Jr \J{^- 



LJo) j 2-K 



/. 



j{uj - UJo) 



27r 



Lo-stable, LSI, 5-continuous 



Lp-stable, LSL 5-continuous 



Causal, LSI 



Output vanishes ai t = to 



/(^) , ^ , N '\ duj 
h /(-aJo)7r()(a; + uJo) e 



~j(^ + ^o) 



27r 



Anti-causal, LSI 



/(o;) - /(-a;o)e-^("+"«)"' ^ do; 



-j(tj + a;o) 



27r 



Lp-stable and decay preserving 



where the On and hm are arbitrary complex coefficients with the normalization constraint ajv = 1. While 

this reminds us of the textbook formula of an ordinary A^th-order differential system, the non-standard 

aspect here is that the driving term is a white noise process w, which is generally not defined pointwise, 

and that we are not imposing any stability constraint. Eq. ( [15] ) thus covers the general case ([8]l where L 

is a shift-invariant operator with the rational transfer function 

w X ^ (j^)^ + aAf-i(ic^)^~^ H h ai(ja^) + Qq ^ -PAf(j^) 

^""^ bMiju;)'' + --- + hijoj)+bo QAiiM' ^ ^ 

The poles of the system, which are the roots of the characteristic polynomial Pn{C) = + aAf-iC"~^ + 

• • • + ao with Laplace variable C G C, are denoted by {an}n=i- While we are not imposing any restriction 

on their locus in the complex plane, we are adopting a special ordering where the purely imaginary roots 



(if present) are coming last. This allows us to factorize the numerator of (16 1 as 

N /N-no \ / \ 



n=l \ n=l / \m=l 
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with aN-no+m = J^m whcrc no is the number of purely-imaginary poles. The operator counterpart of 
this last equation is the decomposition 

Pjv(D) = (D - aild) • • • (D - a^_„Jd) (D - jwild) • • • (D - jojnM) 

^ V ' V ' 

regular part critical part 

which involves a cascade of elementary first-order components. By applying the proper sequence of right- 
inverse operators from Table 1, we can then formally solve the system as in Q. The resulting inverse 
operator is 



L ^ = Iw„o,t„o • • • I<.^i,ti Tlsi 



(18) 



shift- variant 



with 



Tlsi = (D - a^_„Jd)-i • • • (D - aild)-^QM(D), 
which imposes the no boundary conditions 

(D-jw„Jd)s(t)|,^,_^^__^ 

^ {■D-ju;2ld)---{B-jiVnold)s{t)\,^,^ 
The corresponding adjoint operator is given by 








(19) 



L 



-1* 



(20) 



shift-variant 

and is guaranteed to be a continuous linear mapping from S into TZ by Theorem 1, the key point being 
that each of the component operators preserves the rapid decay of the test function to which it is applied. 



The last step is to substitute the explicit form (|20[) of L ^* into (|10[), which yields the characteristic form 



of the stochastic process s defined by ( |T5| ) subject to the boundary conditions ( (T9| ). 

We close this section with a comment about commutativity: while the order of application of the 



operators Qm(D) and (D — a„Id)^^ in the LSI part of (18 1 is immaterial (thanks to the commutativity 
of convolution), it is not so for the inverse operators luj^,to that appear in the "shift- variant" part of the 
decomposition. The latter do not commute and their order of application is tightly linked to the boundary 
conditions. 
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IV. Properties of generalized stochastic processes 

A. Green functions and exponential B-spline calculus 

The foundation of the exponential spline calculus is tiiat we can always factor an A^th-order differential 
operator into a cascade of first-order operators Pq.„ = (D — cunld) where the q;„ (complex poles) are the 
roots of the characteristic polynomial; i.e., 

Piv(D) = + ajv-iD^-i + • • • + aiD + aold 

= -Paw ■ ■ ■ Pqi ~ -f*(ai,...,aiv) 

where the right-hand side concatenated operator notation is self-explanatory. This allows us to express 
the Green function of with pole vector a = (ai, . . . , qat) as the convolution of the Green functions 
of its elementary constituants 

Pa{t) = {Pai * Pa2 ■ ■ ■ * Pa,^){t) (21) 

with 

f u(t)e"* if Re(a) < 

Pait) =1 ^ ^ (22) 

[ -n(-i)e"* else. 

The so-defined Green function pait) is necessarily of slow growth; it specifies the impulse response of 
the LSI inverse operator P~^, which is well-defined over S, 

but not necessarily bounded, as we have seen in the case of purely imaginary poles. 
Next, we observe that by applying the finite difference operator 

Aa/(0 = /W-e"/(t-l) 

to the function Pa(t), we are able to construct a compactly-supported function: the first-order exponential 
B-spline with parameter a 

' l[o,i)(t)e°* if Re(a) <0 
l[o,i)(*)e"^*"'^ else. 

The generaUzation of this scheme yields the iVth-order B-spUne with parameter vector o: = (ai, . . . , a^) 



/3«(i) = A„p«(i) = i/3a^ * /3a. • • • * Pa^){t). (23) 
These functions have the following properties: 
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They are smooth and well-locaUzed: compactly supported in [0, N], bounded, and Holder continuous 
of order — 1. 

They are piecewise-exponential with joining points at the integer and a maximal degree of smoothness 
(spline property). For a = (0, . . . , 0), one recovers Schoenberg's classical polynomial B-splines of 



degree - 1 ||38|, |[39|. 

They are the shortest elementary constituents of splines: the functions {/3a{t — n)}nez forms a Riesz 
basis of the corresponding family of exponential splines with knots at the integers. 



The crucial formula for our purpose is the equivalent operator interpretation of the B -spline formula (23 1: 

AqP^ V = ^aPcx *(p = Pc.*(p, (24) 

which we will now put to good use in order to stationarize and localize the effect of the inverse operators 
that were encovered in the previous sections. 



Theorem 2: Let {I*^ tj\m=i ^^^^ G M be a series of operators of the type defined by ( 13 1 and let 
{A*^ }^'^]^ be some corresponding adjoint localization operators with A*^^(^(t) = (p{t) — e^'^"^Lp{t + \). 
Then, 



for all If ^ S, where /3(jYji,... jaj„(,) exponential B-spline kernel as defined by (23 1. Since the latter is 
bounded and compactly-supported, the resulting convolution operators are BIBO-stable and 5-continuous. 

Proof: First, we observe that A*^^ f{^) = (1 ~ e^'^"'e^'^)f{uj). Using the bottom formula of Table 
I, we then evaluate the Fourier transform of g{t) = I*^^^A*^ /(t) as 



-j{uj + UJm) 



-JUJ - JUJ, 



m 



where we identify the right-hand side factor as /3j(^^ {—'^) where /3a (w) = is the Fourier transform 



of the first-order exponential B-spline with parameter a. This proves that I* ^ A*^ / = /3V * f for 



T* ^ , — 

any LOm , G IR- Using the property that the order of application of stable convolution operators such 
as A*^ can be changed (commutativity), we start with I* ^ A*^^ / and progressively work our way 



outwards to show that I*, . • • • I* . A* , • • • A* , o? = /3V * • • • * * ip, which, thanks to (23 1, 
yields the desired result. The second formula is established in the same way. ■ 
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The interpretation of the second relation is that the difference operators Aj^^^ annihilate the sinusoidal 
components that are in the null space of (D — jtOn) so that the effect of l^i^^t^ becomes indistinguishable 



from that of the non-regularized shift-invariant inverse I^^^^^. By combining this result with (24), we 
obtain a stable LSI substitute for the original inverse operator with the added benefit of a much better 
localization. 

Corollary 1: Let L^^* be the A^th-order (not necessarily shift-invariant) inverse operator specified by 
([20]). Then, 

where /3l is the generalized B -spline kernel 

M 

PUt) = QM(D)/3«(t) = &mD™/3«(t). (25) 

m=l 

The latter is a linear combination of derivatives of the A^th-order exponential B-spline /3a (t) with 
parameter vector a = (ai, . . . , oat), and is therefore compactly-supported over the time-interval [0, N]. 

The intuition behind this result is that we are localizing the system's response by canceling the poles 
of its frequency response; i.e., a pole at juj = a„ is neutralized by a corresponding zero of 1 — e"""-?"^ 
(the frequency response of Aq,^). 

B. Characterization of generalized stochastic processes 

We now proceed with the explicit characterization of a broad class of stochastic processes governed by 
the innovation model ([8]l. Besides the fundamental issue of the solvability of such an operator equation 



which has already been addressed in Section lE-C one needs to make sure that the solutions are bona 



fide generalized stochastic processes. The answer, of course, is dependent upon whether or not we are 
able to exhibit an (adjoint) inverse operator T = L^^* that is sufficiently well-behaved for the resulting 
characteristic form Zs{^) = ^«,(T(^) to be admissible. To that end, we can rely on the following theorem 
whose proof is given in Appendix II. 

Theorem 3 (Admissibility): Let Zs{f) = e^«-^^^^'^ '^^ where the scalar function f{u) is of the generic 
Levy-Khintchine form ^ and where T is an operator acting on 99 G 5. Then, Zs{^) is a continuous, 
positive-definite functional on S such that -Zs(O) = 1, provided that any one of the conditions below is 
met: 

1) T is a continuous linear map from S into itself. 
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2) T is a continuous linear map from S into Lp and /(u) is such that \f{u)\ + \u\ • < C\u\^ 

for all ti G M, where 1 < p < oo and C is a positive constant. 

The simplest scenario is when L^^ is LSI and can be decomposed into a cascade of BIBO-stable and 
ordinary differential operators. If the BIBO-stable part is rapidly-decreasing, then L^^ is guaranteed to 
be 5-continuous. In particular, this covers the case of an A^th-order differential system without any pole 



on the imaginary axis, as justified by our analysis in Section III-C 



Property 1 (Generalized stationary processes): Let L^^ (the right-inverse of some operator L) be a 5- 
continuous convolution operator characterized by its impulse response = L^^S. Then, the generalized 
stochastic processes that are defined by Zs{ip) = exp (/r /(pl * dt) where f{u) is of the generic form 
([3]) are stationary and well-defined solutions of the operator equation ([8]) driven by some corresponding 
white noise process w. 

Proof: The fact that these generalized processes are well-defined is a direct consequence of the 
Minlos-Bochner Theorem since L^^* (the convolution with p^) satisfies the strictest version of admissi- 
bility in Theorem |3] The stationarity property is equivalent to Zs{(p) = Z{ip{- — to)) for all to £ it 
is established by simple change of variable in the inner integral using the basic shift-invariance property 
of convolution; i.e., (p^ * 'fi' — to)) (t) = {Pi^ * ip){t — to). ■ 

The above characterization is not only remarkably concise, but also quite general. It extends the tradi- 
tional theory of stationary Gaussian processes, which corresponds to the choice f{u) = — '^v?. The Gaus- 
sian case results in the simplified form f {L""^* (p) dt = — ^Upl * v\W^ — /r 
(using Parseval's identity) where $5(0;) = is the spectral power density that is associated with 

the innovation model. The interest here is that we get access to a much broader family of non-Gaussian 
processes (e.g., generalized Poisson or alpha-stable) with matched spectral properties since they share 
the same whitening operator L. 

The characteristic form condenses all the statistical information about the process. For instance, by set- 
ting LP = uj6{--to), we can explicitly determine Zs{ip) = {e^^'''^^} = {e^'^^^*")} = T{p{s{to))}{-uj), 
which yields the characteristic function of the first-order probability density of the sample values of the 
process: p(s(to)) = p{s)- In the present stationary scenario, we find thatp(s) = J^~^{exp (Jj^ /( — a;/>L(i)) dt)}(s), 
which requires the evaluation of an integral followed by an inverse Fourier transform. While this type 
of calculation is only tractable analytically in special cases, it may be performed numerically with the 
help of the FFT. Higher-order density functions are accessible as well as at the cost of some multi- 
dimensional inverse Fourier transforms. The same applies for moments which can be obtained through 



a simpler differentiation process, as exemplified in Section IV-C 
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The further reaching aspect of the present formulation is that it is also applicable to the characterization 
of non-stationary processes such as Brownian motion and Levy processes, which are usually treated 
separately from the stationary ones, and that it naturally leads to the identification of a whole variety 
of higher-order extensions. The commonality is that these non-stationary processes can all be derived as 
solutions of an (unstable) A^th-order differential equation with some poles on the imaginary axis. This 



corresponds to the setting in Section III-C with no > 0. 

Property 2 (Generalized Nth-order Levy processes): Let L^^ (the right-inverse of an A^th-order dif- 



ferential operator L) be specified by (18 1 with at least one non-shift-invariant factor \u)i,ti- Then, the 
generalized stochastic processes that are defined by 2s('/') = 6xp (Jj^ /(L~^*(/9) dt), where f{u) is of 
the generic form ^ subject to the constraint + |m| • < C\u\^ for some p > 1, are well- 

defined solutions of the stochastic differential equation ([TSj) driven by some corresponding Levy white 



noise w. These processes satisfy the boundary conditions (19 1 and are non-stationary. 

Note that the required condition on f{u) is satisfied by the great majority of the members of the 
Levy-Kintchine family. For instance in the Poisson case, we can show that \u\ ■ < A|n|£'{|a|} and 

f{u) < A|n|£^{|a|} by using the fact \e^^ — 1| < this implies that the bound in Theorem 3 with p = 1 
is always satisfied provided that the first (absolute) moment of the amplitude PDF Pa{a) is finite. The 
only cases we are aware of that do not fulfill the condition are the alpha-stable noises with < a < 1, 
which are notorious for their exotic behavior. 



Proof: The result is a direct consequence of the analysis in Section III-C — in particular, Eqs. (18l- 
(20l — and Proposition 2. The latter implies that L~^*ip is bounded in all Loo,m norms with m > no 
where no is the number of purely imaginary poles, as specified earlier. Since S C Loo,m C Lp for m > 1 
and the Schwartz topology is the strongest in this chain, we can infer that L~^* is a continuous operator 
from S onto any of the Lp spaces with p > I. The existence claim then follows from the combination of 
Theorem [3] and Minlos-Bochner. Since L^^*ip is not shift-invariant, there is no chance for these processes 
to be stationary, not to mention the fact that they need to fulfill the boundary conditions ( [T9| ). ■ 
Conceptually, we like to view the generalized stochastic processes of Property |2] as "adjusted" versions 
of the stationary ones that include some additional sinusoidal (or polynomial) trends. While the generation 
mechanism of these trends is random, there is a deterministic aspect to it because it imposes the boundary 
conditions ( [T9| ) at ti,--- ,tno- The class of such processes is actually quite rich and the formalism 
surprisingly powerful. We shall illustrate the use of Property 2 in Section V with the simplest possible 
operator L = D which will gets us back to Brownian motion and the celebrated family of Levy processes. 
We shall also show how the well-known properties of Levy processes can be readily deduced from their 
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characteristic form. 

Since the study of Levy processes is usually perceived as presenting a greater level of difficulty than 
the classical theory of stationary processes, one could expect that the same applies to the generalized 
processes of Property 2 because: 1) they are solutions of unstable equations (poles on the imaginary axis), 
and 2) they are non-stationary, which complicates their statistical characterization. While this is indeed 
the case, there is an interesting, practical way to get around the main difficulty by taking advantage 



of the localizing properties of the finite difference operators introduced in Section IV-A These do 
actually provide the proper extension of the classical notion of increments for Brownian motion and Levy 
processes. 

Property 3 (Generalized increment processes): Let s be a generalized stochastic process whose char- 
acteristic form is Zs{<f) = 2^^(L~^*(/?) where Zyj is a white noise functional and where L~^* is given by 



(20 1 (differential system of order N with pole vector a and driving operator Qa/(D) = X]m=i ^mD™). 
The corresponding generalized increment process 

is well-defined and stationary (irrespective of any stability consideration). Its characteristic form is given 



by = * V') where /3l is the generalized B-spline kernel defined by (25 1. 

Proof: : Corollary [T] implies that A^L^^it; = Pi,* w. Since the convolution with the compactly- 
supported kernel /3l defines a continuous LSI operator on S, we can invoke Property 1 with p = Pi,, 
which yields the desired result. ■ 
Since the generalized B-spline /3l is Holder-continuous of order N — M — 1, the above characterization 
allows us to infer that the two processes Sd and s are {N — M — 2) times differentiable in the classical 
sense. In fact, the processes are well-defined pointwise as soon as > M, which is the minimum 
requirement for continuity in the mean-square sense |40|. The other direct implication is that the samples 
of the generalized increment process, Sd{ti) and Sd{t2), are independent as soon as \ti — t2\ > N (due to 
the finite support property of the exponential B-spline /3l)- This means that working with the increment 
process Sd{t) has the remarkable feature of completely suppressing long-range dependencies. 

C. Moments and correlation 

The covariance form of a generalized (complex-valued) process s is defined as: 



I3s{(pi,f2) = ^{{s,ipi) ■ {S,(p2)}. 
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where (s, Lp'i) = {s, 992) when s is real-valued. Thanks to the moment generating properties of the Fourier 
transform, this functional can be calculated from the characteristic form Zs{ip) as 



2 d'^Zs{uJiipi +i02^2) 



(26) 

a;i=0, 0^2=0 



dijjiduj2 

where we are implicitly assuming that the required partial derivative of the characteristic functional exists. 
The autocorrelation of the process is then obtained by making the formal substitution ^pi = 5{- — ti) and 

V?2 = - t2): 

Rs{tl,t2) = <^{s{h)sit2)} = Bs {5{- - h), 6i- - t2)) . 

Alternatively, we can also retrieve the autocorrelation function by invoking the kernel theorem: Bs{y:>i,^2) ■ 

Rs{ti,t2)^l{tl)^p{t2) dti dt2. 

The concept also generalizes for the calculation of the higher-order correlation forrrj^ 

S'{{s,ipi) ■ {s,ip2) ■ ■ ■ {s,(Pn)} = {-jY 



d'^Zs{u}i(pi H \-UJNfN) 



which provides the basis for the determination of higher-order moments and cumulants. 

Here, we concentrate on the calculation of the second-order moments, which happen to be independent 
upon the specific type of noise. For the cases where the covariance is defined and finite, it is not hard to 

is 



show that the generic covariance form of the white noise processes defined in Section II-C 

Bn,{(pi,(p2) = al {lpi,lp2), 

where is a suitable normalization constant that depends on the noise parameters b2,V). We then 
perform the usual adjoint manipulation to transfer the above formula to the filtered version s = L~^w of 
such a noise process. 

Property 4 ( Generalized correlation ): The covariance form of the generalized stochastic process whose 
characteristic form is Zs{(p) = Zw{L~^*ip) where is a white noise functional is given by 

Bsiipi,^2) = ctq {L''^*ipi,L''^*ip2) = ctq (L-iL"^Vi;¥'2), 

and corresponds to the correlation function 

Rsih,t2) = ^{sih) • ^} = {L^L-^*d{- - h),6{- - t2)). 

The latter characterization requires the determination of the impulse response of L^^L^^*. In particular, 
when L^^ is LSI with convolution kernel pL £ -^i> we get that 



*For simplicity, we are only giving the formula for a real-valued process. 
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which confirms that the underlying process is wide-sense stationary. Since the autocorrelation function 
rs(r) is integrable, we also have a one-to-one correspondence with the traditional notion of power 
spectrum: <I>s(a;) = J^{rs}{(jj) = j^^~y|i> where L{u)) is the frequency response of the whitening 
operator L. 

The determination of the correlation function for the non-stationary processes associated with the 



unstable versions of (15 1 is more involved. It can be bypassed if, instead of s{t), we consider the 
generalized increment process Sd{t) = Acts{t) of Property 3. 

Property 5 (Reduction of correlation distances): Let s be a generalized stochastic process whose char- 
acteristic form is Zs{ip) = Zw{L~^*(p) where is a white noise functional and where L~^* is given by 



(20 1 (differential system of order N with pole vector a and driving operator Qm(D) = X]m=i ^mD"*). 



Then, the correlation form of Sd{t) = Aas{t) can be written as 



where /3l is the generahzed B-spline defined by (25 1. The corresponding covariance function is 



RsAh,t2) = S { Aos(*i) • A„s(i2)} = ol * ^l) (t2 - ti) 
which vanishes for (^2 — h) ^ [—N, 

The above result is universal in the sense that it does not distinguish between the stable and unstable 
cases; it can handle A^th-order systems in full generality. 

The conclusion of this section is that one can apply standard techniques from system theory (deter- 
mination of impulse responses) to determine the characteristic form of the whose class of Gaussian and 
non-Gaussian stationary processes defined by ( [T5] ). The functional tools from spline theory (exponential 
B-spline calculus) are helpful as well, especially for the handling and characterization of the non-stationary 
variants of these stochastic processes. 

V. Levy processes revisited 

We now illustrate our method by specifying classical Levy processes — denoted by W{t) — via the 
solution of the (marginally unstable) stochastic differential equation 

-^W{t) = w{t) (27) 

where the driving term w is one of the independent noise processes defined earlier. We shall consider 
the solution W{t) for all t E M, but we shall impose the boundary condition W{to) = with to = to 
make our construction compatible with the classical one which is defined for t > 0. 
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In writing Eq. (27 1, we have intentionally abused the notation in order to be provocative: indeed, when 
w is Gaussian, the solution of this differential equation turns out to be the classical Wiener process, 
which is known to be continuous everywhere (almost surely), but nowhere differentiable in the classical 



sense! Of course, this only makes sense if we interpret (27 1 as {DW, (f) = {w, (p) for all (f E S. 



A. Distributional characterization of Levy processes 



The direct application of the operator formalism developed in Section III yields the solution of (27 1 

rt 



Wit) = lo,owit) 



£w{t)dt, t>0 
-Jl'w{t)dt, t<0 

where Io,o is the unique right inverse of D that imposes the required boundary condition at t = 0. The 
Fourier-based expression of this anti-derivative operator is obtained from the 6th line of Table I by setting 

(wo, to) = (0,0) 

e-''^* - 1 dw 



while its adjoint is specified by 



Io,os(t) = / s(w)- . 



iiMt) = I mzm,!^.^, (28) 

We can then invoke Property 2 to the obtain the characteristic form of the Levy process 

Zwiv) = ^^(lS,o<^) (29) 

which is admissible provided that the Levy function / fulMls the condition in Theorem 3. 

We get the characteristic function of the sample values of the Levy process W{ti) = {W,5{- —ti)), 
by making the substitution ip = u:i5{- — ti) in (29 1: Zvk(iwi(5(- — ti)) = Zw{yJi^*QQ6{- — ti)) with ti > 0. 
Here, we use the Fourier-domain definition (28 1 to evaluate Ioo'^(t ~ *i) = l[o,ii)(*) which happens to 
be an indicator function. Since l[o,ti)(t) is equal to one for t G [0,ti) and zero elsewhere, it is easy to 
evaluate the integral over t in (j2]) where f{u) is given by ([3]), which yields 

= exp{[ f [u;il[o,<o(i)] dt) = exp {t,f [oj^] ) 
Jr 

This result is equivalent to the celebrated Levy-Khinchine representation of the process | [23| , ||4T|. 
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B. Levy increment process 

The Levy increment process is defined as /S.W{t) = W{t) — W{t — 1) where the classical finite- 
difference operator A = A(o) is the discrete counterpart of the whitening operator D. The application 
of Property 3 yields the characteristic form Zu,(/3^^^ * Lp) where 0^Q-^{t) = /3(o)(~*) is the anti-causal 
B-spline of degree 0, which is piecewise-constant and compactly supported in t G (—1,0]. We get the 
pointwise characterization in essentially the same way as before: {e-''^i^^(*i)} = Z/^wiy^i^i' — ^i)) = 
Zw{oJi0^Q-^{- — tif). Here too, we can evaluate the integral f[ujil3^Q^{t — ti)] dt based on the fact that 
the B-spline f^^-^{t — ti) is equal to one for t G (—1 + and zero elsewhere, which leads to the 

particularly simple result 

^|gja.,AH/{*.)} = exp (/(^i)) = pid(u;i), (30) 

where f{u) is given by ([3]). Interestingly, this corresponds to the generic form of the characteristic function 
of a distribution {pid) that is infinitely-divisible. 

To investigate the joint dependencies of the increment process, we turn our attention to the 2-D 
characteristic function of the joint distribution p(AVF(ti), IS.W{t2j) which is obtained by evaluating the 
characteristic form of /SW for ip = uji5{- — ti) +oj25{- — t2): p{u:i,0J2) = '2«)(wi/3^-)(- — ti) +a;2/3^Q-)(- — 
t2)). When the B-splines are non-overlapping, that is, when \ti — t2\ > 1, we take advantage of the 
factorization property of the noise functional to show that 



where pid(wi) is defined in (30). This last factorization result implies that the samples of the increment 
process are independent past the correlation distance of 1 (support of the B-spline), which corresponds 
to the step size of the finite-difference operator. In particular, this implies that the integer samples of 
Al^ are i.i.d. random variables; they correspond to a stationary independent discrete noise process whose 
probability density function pid(AVF) = T^^{pid{—uj)}{AW) may be obtained by taking the inverse 



Fourier transform of (30 1. 



Note, however, that the continuous-domain version of the increment process is not rigorously inde- 
pendent: it will exhibit dependencies over a range corresponding to the support of the B-spline. If the 
driving noise is white with a unit variance, then the autocorrelation function of the increment process is 
given by 

^{AWih) ■ AW{t2)} = (/3^o) * /3(0))(il - ^2) = tTiih - t2) = /3(o,o)(tl - ^2 + 1). 
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poles = [0], zeros = [] 




Fig. 1. Examples of Levy motions W{t) with increasing degrees of sparsity. (a) Brownian motion with Levy triplet (0, 1, 0). (b) 
Levy-Laplace motion with (0,0, ^j^)- (c) Compound Poisson process with (0,0, A^e*" Z^) with A = ^. (d) Symmetric 
Levy flight with (O, 0, l/|al°+^) and a = 1.2. 



C. Examples of Levy processes 

Realizations of four different Levy processes are shown in Fig. 1 together with their Levy triplets 
(61,62,^^(0)) where v{a) is the Levy density such that v{a)da = V{da). The first signal is a Brownian 
motion (a.k.a. Wiener process) that is obtained by integration of a Gaussian white noise. While the sampled 
version of AVF is i.i.d. in all cases, it does not yield a sparse representation in this first instance because 
the underlying distribution remains Gaussian. The second process, which may be termed Levy-Laplace 
motion, is specified by the Levy density v{a) = e~l"l/|a| which is not in Li. By taking the inverse 



Fourier transform of ( |30| ), we can show that its increment process has a Laplace distribution p4| ; note 
that this type of generalized Gaussian model is often used to justify sparsity -promoting signal processing 



techniques based on £1 minimization 1 42 1-| 44 1 . The third piecewise-constant signal is a compound Poisson 
process. It is intrinsically sparse since a good proportion of its increments is zero by construction (with 
probability e~^). The fourth example is an alpha-stable Levy motion (a.k.a. Levy flight) with a = 1.2. 
Here, the distribution of /\W is heavy-tailed (SaS) with unbounded moments for p > a. Although this 
may not be obvious from the picture, this is the sparsest process of the lot because it is ^^-compressible 
in the strongest sense | [37{ . Specifically, we can compress the sequence such as to preserve any prescribed 
portion r < 1 of its average la energy by retaining an arbitrarily small fraction of samples as the length 
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of the signal goes to infinity. 



D. Link with conventional stochastic calculus 



Thanks to (27), we can view a white noise w = W as, the weak derivative of some classical Levy 
processes W{t) which is well-defined pointwise (almost everywhere). This provides us with further 
insights on the range of admissible white noise processes of Section II.C which constitute the driving 
terms of the general stochastic differential equation ([8]). This fundamental observation also makes the 
connection with stochastic calculusp] |[26t, |45|, which avoids the notion of white noise by relying on the 



use of stochastic integrals of the form 

s{t)= I p{t,t')AW{t') 



where is a random (signed) measure associated to some canonical Brownian motion (or, by extension, 
a Levy process) and where p{t, t') is an integration kernel that formally corresponds to our inverse 
operator L~^. 

VL Conclusion 

We have set the foundations of a unifying framework that gives access to the broadest possible class 
of continuous-time stochastic processes that are specifiable by linear, shift-invariant equations, which is 
beneficial for signal processing purposes. To illustrate the method, we have investigated the extended 
family of Levy processes, which, in our view, provide the simplest and most basic examples of sparse 
processes, despite the fact that they are non-stationary (due to their pole at the origin). 

The proposed class of stochastic models and its corresponding mathematical machinery (Fourier 
analysis, characteristic functional, and exponential B-spline calculus) has a large range of applicability, 
as we shall illustrate in the companion paper ||30|. Our formulation has the advantage of maintaining 
backward compatibility with the classical theory of Gaussian stationary processes, while introducing a 
large variety of new processes whose properties are better matched to the currently-dominant paradigm 
in the field which is focused on the notion of sparsity. 

^The Ito integral of conventional stochastic calculus is based on Brownian motion, but the concept can also be generalized 
to Levy driving terms using the more advanced theory of semimartingales |j45J. 
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Appendix I: Positive-definite functionals 



We start by recalling the fundamental notion of positive-definiteness for univariate functions |46|. A 
complex-valued function / of the real variable oj is said to be positive-definite iff. 

N N 
m=l n=l 

for every possible choice of uji, . . . ,ujn € M, ■ ■ ■ ,Cn S C and N G This is equivalent to the 
requirement that the N x N matrix F whose elements are given by [F]mn = /(wm — ^n) is positive 
semi-definite (that is, non-negative definite) for all N, no matter how the a;„'s are chosen. 

Bochner's theorem states that a bounded, continuous function / is positive-definite if and only if it is 
the Fourier transform of a positive and finite Borel measure fj,: 

In particular, Bochner's theorem implies that / is the characteristic function of a random variable — that 
is, f{uj) = S'{e^^'^} = f^e^'^^n{dx) where x is the random variable with probability measure fi — iff. / 
is continuous, positive-definite and /(O) = 1. Note that the above results and formulas also generalize to 
the multivariate setting. 

These concepts carry over as well to functionals on some abstract nuclear space £, the prime example 



being Schwartz's class S of smooth and rapidly-decreasing test functions |19|. 

Definition 1: A complex-valued functional L{ip) defined over the function space £ is said to be 
positive-definite iff. 

N N 

XI XI ^^^^ ~ ^^)^^^n > 
m=l n=l 

for every possible choice of ipi, . . . , Lp^ & £, ^i, . . . , £ C and N G N+. 

Theorem 4 (Minlos-Bochner): Given a functional Z{lp) on a nuclear space £ that is continuous, 
positive-definite and such that Z(0) = 1, there exists a unique probabihty measure on the dual space 
£' such that 



where (s, ip) is the dual pairing map. One further has the guarantee that all finite dimensional probabilities 
densities that can be derived from Z{ip) by setting 99 = uitpi + • • • + ujn^n are mutually compatible. 

The characteristic form therefore uniquely specifies the generalized stochastic process s = s{ip) in 
essentially the same way that the characteristic function fully determines the probability measure of a 
scalar or multivariate random variable. 
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Appendix II: Proof of Theorem 3 



1) As is a generalized random process, Z^^ is a continuous functional on S. This, together with the 
assumption that T is a continuous operator on S, implies that the composed functional Zs{}p) := Zy,{^ip) 
is continuous on S. 

Given the functions (^i, . . . , (/pjv in 5 and some complex coefficients ^i, . . . , Cat, 



l<m,n<N 



l<m,n<N 



l<m,n<N 



>0. 



(by the linearity of the operator T) 
(by the positivity of Zyj) 



This proves the positive-definiteness of the functional Zg on 5. 
Clearly, Z,(0) = Z^{TO) = Z^{0) = 1. 

2) By the assumption on the operator T, Tip G Lp for all ip e S. This together with the assumption 
\f{u)\ < C\u\P impUes that Zs{ip) = e^p{J^f{Tip{t)) dt) is well-defined for all ip e S.By the hnear 
property of the operator T and /(O) = 0, we obtain that Zs{0) = 1. The positive-definiteness of the 
functional Zg is estabUshed by an argument similar to the one used above. Finally we prove the continuity 
of the functional Zg on S: Let {(pn}'^=i be a convergent sequence in S and denote its limit in S by tp. 
Then by the assumption on the linear operator T, Tipn converges to Tcp in Lp; that is, 



lim \\Tipn-Tip\\p = 0. 



(31) 



Next, we observe that 



fit) dt 
tP-^ dt 

<C{\v\P-^ + \u 



<C 



<Cmax(|u|^ "il'wr ^ }\u — v 



V 



p— In 

P-^)\u-v\. 



(by the assumption on /) 



(by the triangle inequality) 



September 1, 2011 



DRAFT 



32 



We then have 

/ f{T^n{t))dt- [ f{T^{t))dt 

Jr Jr 
<C [ \Tipit)r^\Tipnit)-Tipit)\ + \Tipr,{t)-T^it)\Pdt 

JR 

<c(^\\Tip\\P-^\\Tipn-Tip\\p + \\Tipn-Tip\\P^ (by Holder's inequality) 

— >0 as n — ^ oo, (by ( |3T] )) 

which proves the continuity of the functional Zg on 5. ■ 
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